#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <unistd.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <sys/time.h>

#define PHI (24.0f/100.0f)
#define PI (3.14159f)

#define DEG2RAD(d,r) do{r=PI/180*d;}while(0)

#define KNRM "\x1B[0m"
#define KRED "\x1B[31m"
#define KGRN "\x1B[32m"
#define KYEL "\x1B[33m"
#define KBLU "\x1B[34m"
#define KMAG "\x1B[35m"
#define KCYN "\x1B[36m"
#define KWHT "\x1B[37m"

enum times {START, END, DIFF};

void usage(void)
{
	printf("Syntax Error!\n");
	printf("Usage: ./heatdist <iterations> <size> <heat radius> <heat value>\n");
	exit(-1);
}

unsigned long calcRuntime(struct timeval *t)
{
	/* calc time difference */
	t[DIFF].tv_sec = t[END].tv_sec - t[START].tv_sec;
	if (t[START].tv_usec > t[END].tv_usec)
	{
		t[DIFF].tv_usec = (1e6 - t[END].tv_usec) + t[START].tv_usec;
	}
	else
	{
		t[DIFF].tv_usec = t[END].tv_usec - t[START].tv_usec;
	}
	
	return t[DIFF].tv_sec * 1e6 + t[DIFF].tv_usec;
}

int main(int argc, char *argv[])
{
	/* check programm call */
	if (5 != argc) usage();

#if DEBUG == 1
	char color[50];
#endif

	unsigned int i,j,k;
	unsigned long runtime;
	struct timeval time[3];
	
	double deg, rad;
	unsigned int origin;
	int di, dj;

	unsigned int iterations = atoi(argv[1]);
	unsigned int viewsize = atoi(argv[2]);
	unsigned int heatspot = atoi(argv[3]);	// => m
	double heatval = atof(argv[4]);			// => H

	/* sanity check */
	if (iterations < 0 || iterations > 1e6) usage();
	if (viewsize < 0 || viewsize > 1e6) usage();
	if (heatspot < 0 || heatspot >= viewsize) usage();
	if (heatval < 0.0f || heatval > 127.0f) usage();

	/* adjust dimension to fit the desired size */
	unsigned int worksize = viewsize + 4;
	
	/* allocate memory */
	double **x, **xn, **swp, **circle;
	x = (double**) malloc(worksize * sizeof(double*));
	xn = (double**) malloc(worksize * sizeof(double*));
	circle = (double**) malloc(worksize * sizeof(double*));

	for (i=0; i<worksize; i++) 
		x[i] = (double*) malloc(worksize * sizeof(double));
	for (i=0; i<worksize; i++)
		xn[i] = (double*) malloc(worksize * sizeof(double));
	for (i=0; i<worksize; i++) 
		circle[i] = (double*) malloc(worksize * sizeof(double));

	/* init arrays */
	for (i=0; i<worksize; i++) 
		for (j=0; j<worksize; j++)
		{
			x[i][j] = 0.0f;
			xn[i][j] = 0.0f;
			circle[i][j] = 0.0f;
		}

	/* set initial state of heat spot */
	origin = worksize/2;
	for (i=0; i<heatspot; i++)
	{
		/* sweep all radii */
		for (deg = 0; deg <= 360; deg += 0.5f)
		{
			DEG2RAD(deg, rad);
			//printf("%lf %lf\n", deg, rad);
			di = origin + sin(rad)*i;
			dj = origin + cos(rad)*i;
			//printf("%i %i\n", di, dj);
			circle[di][dj] = heatval;
		}
	}

	/* benchmark */
	gettimeofday(&time[START], NULL);
	for (k=0; k<iterations; k++)
	{
		/* keep heat source active */
		for (i=2; i<worksize-2; i++)
		{
			for (j=2; j<worksize-2; j++)
			{
				if (circle[i][j] != 0.0)
					x[i][j] = circle[i][j];
				//if (x[i][j] > 127.0) x[i][j] = 127.0f;
			}
		}

		/* calculate next step */
		for (i=2; i<worksize-2; i++)
		{
			for (j=2; j<worksize-2; j++)
			{
				double *nv = &xn[i][j]; // next value
				double *cv = &x[i][j]; // current value
				double *cvd = &x[i+1][j]; // cv-down 
				double *cvu = &x[i-1][j]; // cv-up
				double *cvr = &x[i][j+1]; // cv-right
				double *cvl = &x[i][j-1]; // cv-left

				*nv = *cv + PHI * (-4.0f * *cv + *cvd + *cvu + *cvr + *cvl);
			}
		}

#if DEBUG == 1
		// uncomment to print output
		for (i=2; i<worksize-2; i++)
		{
			for (j=2; j<worksize-2; j++)
			{
				double *cv = &x[i][j]; // current value

				if (*cv >= 0.0 && *cv < 18.0)
					strcpy(color, KWHT);
				else if (*cv > 18.0 && *cv < 36.0)
					strcpy(color, KCYN);
				else if (*cv > 36.0 && *cv < 54.0)
					strcpy(color, KMAG);
				else if (*cv > 54.0 && *cv < 72.0)
					strcpy(color, KBLU);
				else if (*cv > 72.0 && *cv < 90.0)
					strcpy(color, KGRN);
				else if (*cv > 90.0 && *cv < 108.0)
					strcpy(color, KYEL);
				else if (*cv > 108.0 && *cv <= 127.0)
					strcpy(color, KRED);
				else
					strcpy(color, KNRM);

				//printf("%.1lf ", x[i][j]);
				printf("%s# ", color);
			}
			printf("\n");
		}
		printf("iteration: %i/%i\n", k+1, iterations);
		usleep(1e3*50); // ms
		printf("\n\f");
#endif

		swp = xn;
		xn = x;
		x = swp;
		swp = NULL;
	}
	gettimeofday(&time[END], NULL);

	/* free memory */
	for (i=0; i<worksize; i++)
	{
		free(x[i]);
		free(xn[i]);
		free(circle[i]);
	}
	free(x);
	free(xn);
	free(circle);
	
	runtime = calcRuntime(time);
	printf("%i %i %lf %lf\n",
		viewsize, iterations, (double)runtime/1e6, (double)runtime/1e6/iterations);

	return 0;
}

